This document reports the code for the discovery of COVID-19 subgroups by symptoms and comorbidities, evaluated on the nCov2019 open dataset. This document, as well as the supporting functions required for its execution, are available in the COVID-19-SDE-Tool GitHub repository.
The COVID-19 infectious disease has led since December 2019 to a worldwide pandemic which is still under control measures. Researchers worldwide are making huge efforts aiming to a comprehensive understanding of the COVID-19 and related healthcare treatments. This work shows the preliminary results of a Machine Learning (ML) approach to identify subgroups of COVID-19 patients based on their symptoms and comorbidities, aiming to a better understanding of variability of severity patterns. In this work, we particularly address the variability (or heterogeneity) between distinct sources populating the research repositories, given the potential impact that this variability may have in data science and the generalization of its results.
We analyzed the raw nCov-2019 dataset release at 2020-05-11. The nCoV2019 dataset comprises a collection of publicly available information on worldwide cases confirmed during the ongoing nCoV2019 outbreak. We included those cases were at least one symptom and an outcome were available. Then, we fixed duplicates and homogenized values in outcomes, comorbidities and symptoms. We mapped the latter to ICD-10 terms. The final sample included 170 cases.
Then, we applied a Multiple Correspondence Analysis 3-dimensional embedding of symptoms and outcomes and a hierarchical clustering. The proper number of clusters for both age-independent and age group analyses were selected by supervised inspection of group consistency.
We found clinically meaningful patient subgroups based on symptoms and comorbidities for specific age groups and age-independent analyses. However, the two most prevalent source countries were divided into separate subgroups with different manifestations of severity.
If you use this code please cite:
Carlos Sáez, Nekane Romero, J Alberto Conejero, Juan M García-Gómez. Potential Biases in COVID-19 Machine Learning due to Data Source Variability. Submitted.
If you are interested in collaborating in this work please contact us.
For further exploration of the results please visit our COVID-19 Subgroup Discovery and Exploration Tool (COVID-19 SDE Tool).
Install and load the required packages.
installedPackages = rownames(installed.packages())
if(! "text2vec" %in% installedPackages)
{install.packages("text2vec")}
if(! "plotly" %in% installedPackages)
{install.packages("plotly")}
if(! "factoextra" %in% installedPackages)
{install.packages("factoextra")}
if(! "FactoMineR" %in% installedPackages)
{install.packages("FactoMineR")}
if(! "MVN" %in% installedPackages)
{install.packages("MVN")}
if(! "RColorBrewer" %in% installedPackages)
{install.packages("RColorBrewer")}
if(! "kableExtra" %in% installedPackages)
{install.packages("kableExtra")}
if(! "dplyr" %in% installedPackages)
{install.packages("dplyr")}
if(! "stringr" %in% installedPackages)
{install.packages("stringr")}
library(text2vec)
library(plotly)
library(factoextra)
library(FactoMineR)
library(MVN)
library(RColorBrewer)
library(kableExtra)
library(dplyr)
library(stringr)Source the required functions. These can be found at the COVID-19-SDE-Tool GitHub repository. We recommended to download the whole repository, which includes this document .Rmd file.
Load the nCov2019 dataset at 2020-05-11. Check the variables set if you plan to use other versions of the dataset.
dataDate = '2020-05-11'
untar(paste0('data/latestdata_',dataDate,'.tar.gz'), exdir = "data")
filename = paste0('data/latestdata_',dataDate,'.csv')
data = read.csv2(filename, sep = ",", header = TRUE, na.strings = "", stringsAsFactors = FALSE, dec = '.',
colClasses = c( "character", #id
"character", #age [TO CONVERT TO NUMERIC]
"character", #sex
"character", #city
"character", #province
"character", #country
# "factor", #wuhan.0._not_wuhan.1
"character", #latitude [TO CONVERT TO NUMERIC]
"character", #longitude [TO CONVERT TO NUMERIC]
"character", #geo_resolution
"character", #date_onset_symptoms [TO CONVERT TO DATE]
"character", #date_admission_hospital [TO CONVERT TO DATE]
"character", #date_confirmation [MAIN DATE FOR ANALYSIS (mostly complete), TO CONVERT TO DATE]
"character", #symptoms [TO SPLIT AND PREPROCESS]
"character", #lives_in_wuhan
"character", #travel_history_dates
"character", #travel_history_location
"character", #reported_market_exposure
"character", #additional_information
"character", #chronic_disease_binary
"character", #chronic_disease
"character", #source
"character", #sequence_available
"character", #outcome
"character", #date_death_or_discharge [TO CONVERT TO DATE]
"character", #notes_for_discussion
"character", #location
"character", #admin3
"character", #admin2
"character", #admin1
"character", #country_new
"character", #admin_id
"character", #data_moderator_initials
"character" #travel_history_binary
))Convert and derive variables.
# convert some variable types
data$ageNum = as.numeric(data$age) # Note: Some textual ranges are set as NA (e.g, 16-80), which could be set in another variable
data$ageCat = vector(mode = "character", length = nrow(data))
# make age groups <17, 18-49, 50-64, >65
idsLeq17 = data$ageNum <= 17
ids18to49 = data$ageNum > 17 & data$ageNum <= 49
ids50to64 = data$ageNum > 49 & data$ageNum <= 64
idsGeq65 = data$ageNum > 64
data$ageCat[idsLeq17] = '<17'
data$ageCat[ids18to49] = '18-49'
data$ageCat[ids50to64] = '50-64'
data$ageCat[idsGeq65] = '>65'
# set confirmation date as main date, keep the others as difference to the former
data$date_confirmation = as.Date(data$date_confirmation, "%d.%m.%Y")
data$date_onset_symptoms = as.Date(data$date_onset_symptoms, "%d.%m.%Y")
data$date_admission_hospital = as.Date(data$date_admission_hospital, "%d.%m.%Y")
data$date_death_or_discharge = as.Date(data$date_death_or_discharge, "%d.%m.%Y")
data$date_onset_symptoms_difconf = as.numeric(data$date_onset_symptoms-data$date_confirmation )
data$date_admission_hospital_difconf = as.numeric(data$date_admission_hospital-data$date_confirmation )
data$date_death_or_discharge_difconf = as.numeric(data$date_death_or_discharge-data$date_confirmation )
data$date_death_or_discharge_difadm = as.numeric(data$date_death_or_discharge-data$date_admission_hospital )
data$date_death_or_discharge_difsym = as.numeric(data$date_death_or_discharge-data$date_onset_symptoms )Select subdataset with non-missing symptoms and outcomes (assuming missing chronic diseases mean an absence of them).
Includes the semantic processing of variables through homogenization of clinical terms, and additional data filtering.
outcomeLists = lapply(data_symptoms_outcome$outcome, tolower)
outcomeLists = lapply(outcomeLists, function(x) gsub("\\<discharge\\>|\\<discharged\\>|\\<stable\\>|\\<recovered\\>", "Recovered", x))
outcomeLists = lapply(outcomeLists, function(x) gsub("\\<death\\>|\\<dead\\>|\\<deceased\\>|\\<died\\>", "Deceased", x))
data_symptoms_outcome$outcome2 = sapply(outcomeLists, paste, collapse = " ")
validOutcomes = data_symptoms_outcome$outcome2 %in% c("Recovered","Deceased")
data_symptoms_outcome = data_symptoms_outcome[validOutcomes,]chronicTable = read.csv2('data/chronic_diseases.csv', sep = ";", header = TRUE, na.strings = "", stringsAsFactors = FALSE)
chronicTable = chronicTable[,c(1,4)]
chronicLists = sapply(data_symptoms_outcome$chronic_disease, strsplit, "\\,|;|:|and")
chronicLists = lapply(chronicLists, str_trim)
chronicLists = lapply(chronicLists, tolower)
chronicLists = lapply(chronicLists, str_replace_all, " ", "_")
# write.csv2(unique(unlist(chronicLists)), "chronic_diseases_R.csv", row.names = FALSE, quote = FALSE)
chronicLists = lapply(chronicLists, function(x) chronicTable$group[match(unlist(x), chronicTable$text)])
data_symptoms_outcome$chronic_disease2 = sapply(chronicLists, paste, collapse = " ")# split symptoms texts into vector of (yet unprocessed) symptoms
symptomsLists = sapply(data_symptoms_outcome$symptoms, strsplit, "\\,|;|:|and")
# trim blank spaces
symptomsLists = lapply(symptomsLists, str_trim)
# to lowercase
symptomsLists = lapply(symptomsLists, tolower)
# replace blank spaces with underscore
symptomsLists = lapply(symptomsLists, str_replace_all, " ", "_")
# single replacings of synonyms
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<loss_of_apetite\\>", "anorexia", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<pleuritic_chest_pain\\>", "chest_pain", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<diarrhoea\\>", "diarrhea", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<cardiac_arrythmia\\>", "arrythmia", x))
# multiple replacings of synonyms
# RESPIRATORY
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<acute_rhinitis\\>|\\<acute_pharyngitis\\>|\\<nasal_congestion\\>|\\<nasal_discharge\\>|\\<rhinhorrea\\>|\\<coryza\\>|\\<coriza\\>|\\<runny_nose\\>|\\<running_nose\\>|\\<pharyngeal_dryness\\>|\\<pharyngeal_discomfort\\>|\\<sore_throat\\>|\\<colds\\>|\\<dysphagia\\>|\\<dry_throat\\>|\\<throat_discomfort\\>|\\<pharyngalgia\\>|\\<flu_like_symptoms\\>|\\<cold\\>|\\<colds\\>", "acute_nasopharyngitis", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<coughing|.+cough", "cough", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<acute_respiratory_distress\\>|\\<acute_respiratory_distress_syndrome\\>|<acute_respiratory_disease_syndrome\\>", "ards", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<difficulty_breathing\\>|\\<shortness_of_breath\\>|\\<chest_distress\\>|\\<chest_myalgia\\>|\\<chest_tightness\\>|\\<chest_discomfort\\>\\<chest_fatigue\\>|\\<respiratory_problems\\>|\\<gasp\\>|\\<grasp\\>|\\<breathing_difficulty\\>|\\<respiratory_complaints\\>", "dyspnea", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<sputum\\>|\\<little_sputum\\>|\\<phlegm\\>|\\<little_expectoration\\>", "expectoration", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<severe_acute_respiratory_infection\\>|\\<severepneumonia\\>", "pneumonia severe_pneumonia", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<pneumonitis\\>|\\<mild_pneumonia\\>|\\<lesions_on_chest_radiographs\\>", "pneumonia", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<acute_respiratory_failure\\>|\\<hypoxia\\>|\\<hypercapnia\\>", "respiratory_failure", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<acute_respiratory_disease\\>", "unspecified_respiratory_disease", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<sepsis\\>", "\\<septic_shock\\>", x))
# NON RESPIRATORY
symptomsLists = lapply(symptomsLists, function(x) gsub("fever.+|\\<chills\\>|.+chill|\\<sensation_of_chill\\>", "fever", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<acute_renal_failure\\>|\\<prerenal_failure\\>|\\<kidney_failure\\>", "acute_kidney_injury", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<eye_irritation\\>|\\<red_eye\\>", "conjunctivitis", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<joint_pain\\>|\\<joint_tenderness\\>", "arthralgia", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<dizziness\\>|\\<transientfatigue\\>|\\<transient_fatigue\\>|\\<discomfort\\>|\\<body_malaise\\>|\\<exhaustion\\>|\\<feeling_ill\\>|\\<general_malaise\\>|\\<general_weakness\\>|\\<lack_of_energy\\>|\\<lethargy\\>|\\<malaise\\>|\\<systemic_weakness\\>|\\<iredness\\>|\\<general_weakness\\>|\\<weak\\>|\\<weakness\\>|\\<weaknessness\\>|\\<fatigure\\>", "fatigue", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("+.myalgia|\\<myalgias\\>|\\<aching_muscles\\>|\\<backache\\>|\\<bone_pain\\>|\\<body_ache\\>|\\<milagia\\>|\\<mialgia\\>|\\<muscle_soreness\\>|\\<muscular_soreness\\>|\\<muscle_pain\\>|\\<musculoskeletal_pain\\>|\\<muscular_stiffness\\>|\\<pain\\>", "myalgia", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<nausea\\>|\\<vomiting\\>|\\<vomits\\>|\\<emesis\\>", "nausea_vomiting", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<drowsiness\\>|\\<somnolence\\>|\\<obnubilation\\>", "altered_conciousness_mild", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<congestive_heart_failure\\>|\\<myocardial_dysfunction\\>", "heart_failure", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<infarction\\>|\\<myocardial_infarction\\>|\\<acute_myocardial_infarction\\>", "acute_coronary_syndrome", x))
# symptomsLists = lapply(symptomsLists, function(x) gsub("\\<asymptomatic\\>|\\<none\\>", "NA", x))
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<asymptomatic\\>|\\<none\\>|\\<afebrile\\>", "", x))
# deletion of not informative texts
symptomsLists = lapply(symptomsLists, function(x) gsub("\\<19_related_symptoms\\>|\\<covid\\>|\\<severe\\>|\\<between_others\\>|\\<mild\\>|\\<moderate\\>|\\<hypertension\\>", "", x))
#Note: further use of modulators could be revised
#Note: patients with severe only have that value
symptomsLists = lapply(symptomsLists, function(x) x[sapply(x, str_length)>0])
# join again in a new column for automatic word vector creation
data_symptoms_outcome$symptoms2 = sapply(symptomsLists, paste, collapse = " ")Filter by valid age and remove duplicated entries.
To continue with the remaining code please select one of the following chunks of code to execute depending on which Age group you want to analyze. These chunks include the selection of a subset of data for the corresponding age, and setting k as the best number of clusters for each group after a supervised expert review. For the following results we ran the “All ages” group.
Analysis with all ages:
Analysis with ages >65 years:
Analysis with ages between 50-64 years:
Analysis with ages between 18-49 years:
To perform the tokenization of symptoms and comorbidities, provided as a list, we use the performWord2Vec function of the repository.
# Word vectors with word2vec for symptoms
resultsWord2VecSymptoms = performWord2Vec(dataExperiment, targetVariable = 'symptoms2')
# filter out of vocabulary individuals
symptomsVector = resultsWord2VecSymptoms$textsVector[resultsWord2VecSymptoms$idsValidTexts,]
dataExperiment = dataExperiment[resultsWord2VecSymptoms$idsValidTexts,]
## print those symptoms not included to improve filtering (as a check)
# data_symptoms_outcome[!resultsWord2VecSymptoms$idsValidTexts,"symptoms"]
# Word vectors with word2vec for chronic diseases
resultsWord2VecChronic = performWord2Vec(dataExperiment, targetVariable = 'chronic_disease2')
# filter out of vocabulary individuals
chronicsVector = resultsWord2VecChronic$textsVector[resultsWord2VecChronic$idsValidTexts,]
dataExperiment = dataExperiment[resultsWord2VecChronic$idsValidTexts,]
symptomsVector = symptomsVector[resultsWord2VecChronic$idsValidTexts,]Prepare data for analysis, joining symptoms and comorbidities cin a single data.frame and selecting a set of metadata to complement results.
symptomsVector = data.frame(symptomsVector)
colnames(symptomsVector) = paste0("S_",colnames(symptomsVector))
chronicsVector = data.frame(chronicsVector)
colnames(chronicsVector) = paste0("C_",colnames(chronicsVector))
data_analysis = cbind(symptomsVector, chronicsVector)
data_analysis = sapply(data_analysis, as.logical)
data_analysis_metadata = dataExperiment[,c("ageNum", "ageCat", "sex", "country", "source", "date_death_or_discharge_difsym")]
data_analysis_metadata$Outcome = dataExperiment$outcome2Perform Multiple Correspondence Analysis on 3 dimensions.
Perform hierarchical clustering and cut the resultant tree using the best number of clusters k defined above.
distEuclideanMCA = dist(ind$coord)
clusterdistEuclideanMCA = hclust(distEuclideanMCA, method = "ward.D2")
groups <- cutree(clusterdistEuclideanMCA, k=bestk)
sizes <- dataExperiment$ageNum
resultsClustering = list("ind" = ind, "clusterdistEuclideanMCA" = clusterdistEuclideanMCA, "k" = bestk, "groups" = groups, "sizes" = sizes)In the following results the information about the subgroups is consistent within scatter plots, histograms and detailed table, and correspond to the subgroups found in the previous application of clustering to the selected Age group with the corresponding number of clusters k.
Perform first some preparation to facilitate visualizations.
# selection of true symptoms and comorbidities values from MCA
trues = endsWith(names(res.mca$call$marge.col), "TRUE")
coord_T = var$coord[trues,]
rownames(coord_T) = substr(rownames(coord_T),1,nchar(rownames(coord_T))-5)
coord_T_S = coord_T[1:ncol(symptomsVector),]
rownames(coord_T_S) = substr(rownames(coord_T_S),3,nchar(rownames(coord_T_S)))
coord_T_C = coord_T[-(1:ncol(symptomsVector)),]
rownames(coord_T_C) = substr(rownames(coord_T_C),3,nchar(rownames(coord_T_C)))
# configuration of scatter properties
markerSize2d = 35
markerSize2ddiamond = 25
markerSize3d = 300
markerSize3ddiamond = 250
axx <- list(
title = "1<sup>st</sup> component"
)
axy <- list(
title = "2<sup>nd</sup> component"
)
axz <- list(
title = "3<sup>rd</sup> component"
)
# color resources
colorsVars = brewer.pal(n = 2, name = "BrBG")pClusters2d <- plot_ly(x = resultsClustering$ind$coord[,1], y = resultsClustering$ind$coord[,2],
color = as.factor(paste("Subgroup",format(resultsClustering$groups, digits=2))),
size = I(markerSize2d), type = "scatter", mode = "markers",
text = paste0("Patient ID: ",rownames(resultsClustering$ind$coord),"\n Age:",resultsClustering$sizes),
marker = list(sizemode = 'area'),
scene = 'sceneClusters') %>%
layout(title = "Subgroups (obtained at 3D, switch to 3D for better display)", xaxis = axx, yaxis = axy) %>%
config(displaylogo = FALSE)
pClusters2dpClusters3d <- plot_ly(x = resultsClustering$ind$coord[,1], y = resultsClustering$ind$coord[,2], z = resultsClustering$ind$coord[,3],
color = as.factor(paste("Subgroup",format(resultsClustering$groups, digits=2))),
size = I(markerSize3d), type = "scatter3d", mode = "markers",
text = paste0("Patient ID: ",rownames(resultsClustering$ind$coord),"\n Age:",resultsClustering$sizes),
marker = list(sizemode = 'area'),
scene = 'sceneClusters') %>%
layout(title = "Subgroups", scene = list(xaxis=axx,yaxis=axy,zaxis=axz)) %>%
config(displaylogo = FALSE)
pClusters3dpOutcome2d <- plot_ly(x = resultsClustering$ind$coord[,1], y = resultsClustering$ind$coord[,2],
color = as.factor(data_analysis_metadata$Outcome),
# colors = brewer.pal(n = 2, name = "Dark2"),
colors = c("#99004C", "#004C99"),
size = I(markerSize2ddiamond),
text = paste0("Patient ID: ",rownames(resultsClustering$ind$coord),"\n Age:",resultsClustering$sizes),
marker = list(sizemode = 'area', symbol = "diamond"),
scene = 'sceneOutcome') %>%
layout(title = "Outcome", xaxis = axx, yaxis = axy) %>%
add_markers() %>%
config(displaylogo = FALSE)
pOutcome2dpOutcome3d <- plot_ly(x = resultsClustering$ind$coord[,1], y = resultsClustering$ind$coord[,2], z = resultsClustering$ind$coord[,3],
color = as.factor(data_analysis_metadata$Outcome),
colors = c("#99004C", "#004C99"),
type = "scatter3d", mode = "markers",
# size = I(resultsClustering$sizes*20),
size = I(markerSize3ddiamond),
text = paste0("Patient ID: ",rownames(resultsClustering$ind$coord),"\n Age:",resultsClustering$sizes),
marker = list(sizemode = 'area', symbol = "diamond"),
scene = 'sceneOutcome') %>%
layout(title = "Outcome", scene = list(xaxis=axx,yaxis=axy,zaxis=axz)) %>%
config(displaylogo = FALSE)
pOutcome3dpCountry2d <- plot_ly(x = resultsClustering$ind$coord[,1], y = resultsClustering$ind$coord[,2],
color = as.factor(data_analysis_metadata$country),
colors = brewer.pal(n = length(unique(data_analysis_metadata$country)), name = "Set1"),
size = I(markerSize2d), type = "scatter", mode = "markers",
text = paste0("Patient ID: ",rownames(resultsClustering$ind$coord),"\n Age:",resultsClustering$sizes),
marker = list(sizemode = 'area'),
scene = 'sceneCountry') %>%
layout(title = "Country", xaxis = axx, yaxis = axy) %>%
config(displaylogo = FALSE)
pCountry2dpCountry3d <- plot_ly(x = resultsClustering$ind$coord[,1], y = resultsClustering$ind$coord[,2], z = resultsClustering$ind$coord[,3],
color = as.factor(data_analysis_metadata$country),
colors = brewer.pal(n = length(unique(data_analysis_metadata$country)), name = "Set1"),
type = "scatter3d", mode = "markers",
# size = I(resultsClustering$sizes*20),
size = I(markerSize3d),
text = paste0("Patient ID: ",rownames(resultsClustering$ind$coord),"\n Age:",resultsClustering$sizes),
marker = list(sizemode = 'area'),
scene = 'sceneCountry') %>%
layout(title = "Country", scene = list(xaxis=axx,yaxis=axy,zaxis=axz)) %>%
config(displaylogo = FALSE)
pCountry3dCalculate required statistics.
alphaci = 0.05
uniqueGroups = unique(resultsClustering$groups)
nSubgroups = length(uniqueGroups)
resultsBySubgroup = vector("list",length(uniqueGroups))
# colnames(data_analysis) = substr(colnames(data_analysis),3,nchar(colnames(data_analysis)))
for (i in 1:length(uniqueGroups)){
patientGroupIdx = resultsClustering$groups %in% uniqueGroups[i]
nPatientsGroup = sum(patientGroupIdx)
data_analysis_subgroup = data_analysis[patientGroupIdx,,drop = FALSE]
nind = nrow(data_analysis_subgroup)
resultsSymptoms = sapply(data.frame(data_analysis_subgroup[,1:ncol(symptomsVector), drop = FALSE]), function(x) pCI(x, alphaci)*100)
colnames(resultsSymptoms) = str_replace_all(colnames(resultsSymptoms),"\\.", " ")
resultsSymptomsErr = resultsSymptoms[3,] - resultsSymptoms[1,]
resultsComorbidities = sapply(data.frame(data_analysis_subgroup[,-(1:ncol(symptomsVector)), drop = FALSE]), function(x) pCI(x, alphaci)*100)
colnames(resultsComorbidities) = str_replace_all(colnames(resultsComorbidities),"\\.", " ")
resultsComorbiditiesErr = resultsComorbidities[3,] - resultsComorbidities[1,]
# sex, age, recovered statistics
data_analysis_metadata_subgroup = data_analysis_metadata[patientGroupIdx,, drop = FALSE]
femaleStats = pCI(as.character(data_analysis_metadata_subgroup$sex) %in% 'female',alphaci)*100
ageMean = mean(data_analysis_metadata_subgroup$ageNum)
ageErr = qnorm(0.975)*sd(data_analysis_metadata_subgroup$ageNum)/sqrt(nPatientsGroup)
ageStats = c(ageMean, ageErr)
recoveredStats = pCI(data_analysis_metadata_subgroup$Outcome == 'Recovered',alphaci)*100
subgroupResults = list(name = paste("Subgroup", i), symptoms = resultsSymptoms, comorbidities = resultsComorbidities, symptomsErr = resultsSymptomsErr, comorbiditiesErr = resultsComorbiditiesErr, nPatientsGroup = nPatientsGroup, femaleStats = femaleStats, recoveredStats = recoveredStats, ageStats = ageStats)
resultsBySubgroup[[i]] <- subgroupResults
}resultsStats = lapply(resultsBySubgroup, function(x) x$symptoms)
resultsStatsMean = unlist(lapply(resultsStats, function(x) x[1,]))
resultsStatsErr = unlist(lapply(resultsBySubgroup, function(x) x$symptomsErr))
subgroupIds = 1:length(resultsBySubgroup)
groupColors = colorRampPalette(brewer.pal(nSubgroups,"Set2"))(nSubgroups)
i = 1
pSymptomsBar <- plot_ly(x = colnames(resultsStats[[i]]), y = resultsStats[[i]][1,], type = 'bar', name = paste("Subgroup", i),
marker = list(color = groupColors[i]),
error_y = list(type = "data",
array = pmin(100-resultsStats[[i]][1,],resultsBySubgroup[[i]]$symptomsErr),
arrayminus = pmin(resultsStats[[i]][1,],resultsBySubgroup[[i]]$symptomsErr),
color = '#888888'))
for (i in 2:length(uniqueGroups)){
pSymptomsBar <- pSymptomsBar %>% add_trace(x = colnames(resultsStats[[i]]), y = resultsStats[[i]][1,], type = 'bar', name = paste("Subgroup", i),
marker = list(color = groupColors[i]),
error_y = list(type = "data",
array = pmin(100-resultsStats[[i]][1,],resultsBySubgroup[[i]]$symptomsErr),
arrayminus = pmin(resultsStats[[i]][1,],resultsBySubgroup[[i]]$symptomsErr),
color = '#888888'))
}
# config(displayModeBar = FALSE) %>%
pSymptomsBar <- pSymptomsBar %>% layout(title = sprintf("Symptoms by subgroup"),
xaxis = list(title = "Symptoms"),
yaxis = list(title = "% (CI 95%)"),
barmode = 'group'
) %>%
config(displaylogo = FALSE)
pSymptomsBarresultsStats = lapply(resultsBySubgroup, function(x) x$comorbidities)
resultsStatsMean = unlist(lapply(resultsStats, function(x) x[1,]))
resultsStatsErr = unlist(lapply(resultsBySubgroup, function(x) x$comorbiditiesErr))
subgroupIds = 1:length(resultsBySubgroup)
groupColors = colorRampPalette(brewer.pal(nSubgroups,"Set2"))(nSubgroups)
i = 1
pComorbiditiesBar <- plot_ly(x = colnames(resultsStats[[i]]), y = resultsStats[[i]][1,], type = 'bar', name = paste("Subgroup", i),
marker = list(color = groupColors[i]),
error_y = list(symmetric = FALSE,
type = "data",
array = pmin(100-resultsStats[[i]][1,],resultsBySubgroup[[i]]$comorbiditiesErr),
arrayminus = pmin(resultsStats[[i]][1,],resultsBySubgroup[[i]]$comorbiditiesErr),
color = '#888888'))
for (i in 2:length(uniqueGroups)){
pComorbiditiesBar <- pComorbiditiesBar %>% add_trace(x = colnames(resultsStats[[i]]), y = resultsStats[[i]][1,], type = 'bar', name = paste("Subgroup", i),
marker = list(color = groupColors[i]),
error_y = list(type = "data",
array = pmin(100-resultsStats[[i]][1,],resultsBySubgroup[[i]]$comorbiditiesErr),
arrayminus = pmin(resultsStats[[i]][1,],resultsBySubgroup[[i]]$comorbiditiesErr),
color = '#888888'))
}
# config(displayModeBar = FALSE) %>%
pComorbiditiesBar <- pComorbiditiesBar %>% layout(title = sprintf("Comorbidities by subgroup"),
xaxis = list(title = "Comorbidities"),
yaxis = list(title = "% (CI 95%)"),
barmode = 'group'
) %>%
config(displaylogo = FALSE)
pComorbiditiesBarpAges <- plot_ly(x = resultsClustering$groups, y = data_analysis_metadata$ageNum, split = paste("Subgroup", resultsClustering$groups),
type = 'box',
color = groupColors
)
pAges <- plot_ly(x = resultsClustering$groups, y = data_analysis_metadata$ageNum,
type = 'box',
color = paste("Subgroup", resultsClustering$groups),
colors = groupColors
)
pAges <- pAges %>%
layout(
title = "Age by subgroup",
xaxis = list(
title = "Subgroup"
),
yaxis = list(
title = "Age",
zeroline = F
)
) %>%
config(displaylogo = FALSE)alphaci = 0.05
uniqueGroups = unique(resultsClustering$groups)
resultsTable = data.frame(matrix(nrow = ncol(data_analysis)+6, ncol = length(uniqueGroups)))
colnames(data_analysis) = substr(colnames(data_analysis),3,nchar(colnames(data_analysis)))
colnames(data_analysis) = make.names(colnames(data_analysis), unique = TRUE)
colnames(data_analysis) = str_replace_all(colnames(data_analysis),"_", " ")
# Note: this (whihout make.names) returns an error when same text is in symptoms and comorbidities, a solution might be using a column for names, or avoid repated terms
rownames(resultsTable) <- c(sprintf('No. of individuals (n<sub>total</sub> = %d)', nrow(data_analysis)), colnames(data_analysis), 'Females','Age','Recovered','% valid survival (deceased)','Survival days (deceased)')
colnames(resultsTable) <- paste('Subgroup',uniqueGroups)
for (i in 1:length(uniqueGroups)){
patientGroupIdx = resultsClustering$groups %in% uniqueGroups[i]
nPatientsGroup = sum(patientGroupIdx)
# comorbidity statistics
data_analysis_subgroup = data_analysis[patientGroupIdx,,drop = FALSE]
nind = nrow(data_analysis_subgroup)
data_analysis_subgroupT = t(data_analysis_subgroup)
resultsS = sapply(data.frame(data_analysis_subgroup), function(x) pCI(x, alphaci)*100)
subgroupResultColumn = apply(resultsS, 2, function(x) sprintf('%.2f (%.2f-%.2f)',x[1],x[2],x[3]))
# sex, age, recovered statistics
data_analysis_metadata_subgroup = data_analysis_metadata[patientGroupIdx,]
#sexProportions = table(corpusSex[patientGroupIdx])/nPatientsGroup
femaleStats = pCI(as.character(data_analysis_metadata_subgroup$sex) %in% 'female',alphaci)*100
# ageMean = mean(data_analysis_metadata_subgroup$ageNum)
ageStats = mCI(data_analysis_metadata_subgroup$ageNum, alphaci)
recoveredStats = pCI(data_analysis_metadata_subgroup$Outcome == 'Recovered',alphaci)*100
survivalStats = mCI(data_analysis_metadata_subgroup$date_death_or_discharge_difsym[data_analysis_metadata_subgroup$Outcome == 'Deceased'], alphaci, na.rm = TRUE)
# compile final result column
subgroupResultColumn = c(as.character(nPatientsGroup), subgroupResultColumn, sprintf('%.2f (%.2f-%.2f)',femaleStats[1],femaleStats[2],femaleStats[3]), sprintf('%.2f (%.2f-%.2f)',ageStats[1],ageStats[2],ageStats[3]), sprintf('%.2f (%.2f-%.2f)',recoveredStats[1],recoveredStats[2],recoveredStats[3]))
#names(subgroupResultColumn)[c(1,length(subgroupResultColumn)-1,length(subgroupResultColumn))] <- c('n','Females','Median birth date')
# subgroupResultColumn = c(subgroupResultColumn, sprintf('%.2f',sum(!is.na(data_analysis_metadata_subgroup$date_death_or_discharge_difsym))/nind), sprintf('%.2f', mean(data_analysis_metadata_subgroup$date_death_or_discharge_difsym, na.rm = TRUE)))
subgroupResultColumn = c(subgroupResultColumn,
sprintf('%.2f',100*sum(!is.na(data_analysis_metadata_subgroup$date_death_or_discharge_difsym[data_analysis_metadata_subgroup$Outcome == 'Deceased']))/sum(data_analysis_metadata_subgroup$Outcome == 'Deceased')),
# sprintf('%.2f', mean(data_analysis_metadata_subgroup$date_death_or_discharge_difsym[data_analysis_metadata_subgroup$Outcome == 'Deceased'], na.rm = TRUE)))
sprintf('%.2f (%.2f-%.2f)',survivalStats[1],survivalStats[2],survivalStats[3]))
resultsTable[i] <- subgroupResultColumn
}
# names(resultsTable) <- cell_spec(names(resultsTable), background = "yellow")
tTable = kable(resultsTable, format = "html", escape = F) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F) %>%
pack_rows("Symptoms (%, CI 95%)", 2, 1+ncol(symptomsVector)) %>%
pack_rows("Comorbidities (%, CI 95%)", 2+ncol(symptomsVector), 2+ncol(symptomsVector)+ncol(chronicsVector)-1) %>%
pack_rows("Demographics (%|x, CI 95%)", 2+ncol(symptomsVector)+ncol(chronicsVector), 2+ncol(symptomsVector)+ncol(chronicsVector)+1) %>%
pack_rows("Outcomes (%|x, CI 95%)", 2+ncol(symptomsVector)+ncol(chronicsVector)+2, 2+ncol(symptomsVector)+ncol(chronicsVector)+3+1)
groupColors = brewer.pal(n = ncol(resultsTable), name = "Set2")
for (i in 1:ncol(resultsTable)) {
tTable = tTable %>% column_spec(i+1, background = groupColors[i], include_thead = TRUE) %>%
column_spec(i+1, background = "inherit")
}
tTable| Subgroup 1 | Subgroup 2 | Subgroup 3 | Subgroup 4 | Subgroup 5 | Subgroup 6 | |
|---|---|---|---|---|---|---|
| No. of individuals (ntotal = 159) | 54 | 53 | 10 | 23 | 10 | 9 |
| Symptoms (%, CI 95%) | ||||||
| headache | 7.41 (0.42-14.39) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 4.35 (-3.99-12.68) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| severe pneumonia | 0.00 (0.00-0.00) | 3.77 (-1.36-8.90) | 20.00 (-4.79-44.79) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 11.11 (-9.42-31.64) |
| acute coronary syndrome | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 66.67 (35.87-97.46) |
| myalgia | 9.26 (1.53-16.99) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 4.35 (-3.99-12.68) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| acute kidney injury | 0.00 (0.00-0.00) | 3.77 (-1.36-8.90) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 50.00 (19.01-80.99) | 0.00 (0.00-0.00) |
| heart failure | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 20.00 (-4.79-44.79) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 55.56 (23.09-88.02) |
| dyspnea | 5.56 (-0.55-11.67) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 34.78 (15.32-54.25) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| fatigue | 16.67 (6.73-26.61) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 26.09 (8.14-44.03) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| acute nasopharyngitis | 22.22 (11.13-33.31) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 8.70 (-2.82-20.21) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| septic shock | 0.00 (0.00-0.00) | 16.98 (6.87-27.09) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 100.00 (100.00-100.00) | 0.00 (0.00-0.00) |
| respiratory failure | 0.00 (0.00-0.00) | 33.96 (21.21-46.71) | 30.00 (1.60-58.40) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 22.22 (-4.94-49.38) |
| ards | 0.00 (0.00-0.00) | 35.85 (22.94-48.76) | 60.00 (29.64-90.36) | 4.35 (-3.99-12.68) | 70.00 (41.60-98.40) | 11.11 (-9.42-31.64) |
| cough | 59.26 (46.15-72.36) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 56.52 (36.26-76.78) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| pneumonia | 1.85 (-1.74-5.45) | 83.02 (72.91-93.13) | 60.00 (29.64-90.36) | 0.00 (0.00-0.00) | 60.00 (29.64-90.36) | 88.89 (68.36-109.42) |
| fever | 81.48 (71.12-91.84) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) | 78.26 (61.40-95.12) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| Comorbidities (%, CI 95%) | ||||||
| other | 0.00 (0.00-0.00) | 3.77 (-1.36-8.90) | 0.00 (0.00-0.00) | 13.04 (-0.72-26.81) | 0.00 (0.00-0.00) | 11.11 (-9.42-31.64) |
| respiratory | 0.00 (0.00-0.00) | 5.66 (-0.56-11.88) | 0.00 (0.00-0.00) | 26.09 (8.14-44.03) | 40.00 (9.64-70.36) | 0.00 (0.00-0.00) |
| renal | 0.00 (0.00-0.00) | 22.64 (11.37-33.91) | 0.00 (0.00-0.00) | 4.35 (-3.99-12.68) | 10.00 (-8.59-28.59) | 33.33 (2.54-64.13) |
| cardiovascular | 0.00 (0.00-0.00) | 16.98 (6.87-27.09) | 0.00 (0.00-0.00) | 26.09 (8.14-44.03) | 10.00 (-8.59-28.59) | 0.00 (0.00-0.00) |
| metabolic | 0.00 (0.00-0.00) | 58.49 (45.22-71.76) | 10.00 (-8.59-28.59) | 47.83 (27.41-68.24) | 40.00 (9.64-70.36) | 44.44 (11.98-76.91) |
| hypertension | 0.00 (0.00-0.00) | 90.57 (82.70-98.44) | 0.00 (0.00-0.00) | 65.22 (45.75-84.68) | 20.00 (-4.79-44.79) | 55.56 (23.09-88.02) |
| na | 100.00 (100.00-100.00) | 5.66 (-0.56-11.88) | 90.00 (71.41-108.59) | 4.35 (-3.99-12.68) | 60.00 (29.64-90.36) | 11.11 (-9.42-31.64) |
| Demographics (%|x, CI 95%) | ||||||
| Females | 33.33 (20.76-45.91) | 26.42 (14.55-38.28) | 30.00 (1.60-58.40) | 30.43 (11.63-49.24) | 30.00 (1.60-58.40) | 22.22 (-4.94-49.38) |
| Age | 48.19 (43.29-53.09) | 68.66 (65.64-71.68) | 56.50 (47.83-65.17) | 71.65 (66.96-76.34) | 69.70 (62.08-77.32) | 71.56 (62.38-80.73) |
| Outcomes (%|x, CI 95%) | ||||||
| Recovered | 68.52 (56.13-80.91) | 0.00 (0.00-0.00) | 10.00 (-8.59-28.59) | 8.70 (-2.82-20.21) | 0.00 (0.00-0.00) | 0.00 (0.00-0.00) |
| % valid survival (deceased) | 76.47 | 9.43 | 22.22 | 100.00 | 10.00 | 0.00 |
| Survival days (deceased) | 16.00 (11.40-20.60) | 14.40 (11.66-17.14) | 16.00 (12.08-19.92) | 14.33 (11.27-17.39) | 14.00 (NA-NA) | NaN (NaN-NaN) |